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ABSTRACT 

We present a method for simulating numerically the effect of the adiabatic growth 
of black holes on the structure of elliptical galaxies. Using a parallel self-consistent 
held code, we add black holes to N-body realizations of model distribution functions 
for spherical galaxies, using a continuous mass-spectrum. The variable particle mass, 
combined with a simple multiple timestep integration scheme, makes it possible to evolve 
the models for many dynamical times with N ~ fO 6 — fO 8 , allowing high spatial and 
mass resolution. This paper discusses verification of the code using analytic models for 
spherical galaxies, comparing our numerical results of the effect of central black holes 
on the structure of the galaxies with previously published models. The intrinsic and 
projected properties of the final particle distribution, including higher order moments 
of the velocity distribution, permit comparison with observed characteristics of real 
galaxies, and constrain the masses of any central black holes present in those galaxies. 
Our technique is promising and is easily extended to axisymmetric and triaxial galaxies. 
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1. Introduction 

Theories of energy production from quasars and 
active galactic nuclei predict that many present-day 
galaxies contain central massive black holes (MBHs) 
with masses Mbh £ 1O 7 M (Rees 1990). St rong 
observational efforts have revealed many candidate 
galaxies suspected of harboring MBHs, but conclu- 
sive proof remains elusive (see reviews by Dressier 

1989, Gerhard 1992, Kormendy 1992). In tandem 
with past observational studies, and proposed stud- 
ies to be made with the refurbished Hubble Space 
Telescope, Keck and other ground based observato- 
ries (eg. Sargent et al., 1978, Dressier k Richstone 

1990, Lauer et al., 1992a, Lauer et al., 1992b, Sti- 
avelli et al., 1993, Crane et al., 1993, Harms et al., 
1994, van der Marel et al., 1994, Kormendy et al., 
1994), there has been intense theoretical effort to con- 
strain the masses of these claimed central black holes 
and provide theoretical predictions of their conse- 
quences (eg. Bahcall k Wolf 1976, Young 1980, Dun- 
can k Wheeler 1980, Binney k Mamon 1982, Tonry 
1983, Richstone k Tremaine 1985, Shapiro 1985, Bin- 
ney k Petit 1989, Lee k Goodman 1989). 

Using N-body simulations and analytic techniques, 
it is possible to investigate the effects of a central 
MBH on the dynamics of stars in galaxies, and predict 
the range of observational properties of real galaxies 
containing MBHs. A variety of approaches have been 
used to model galaxies containing MBHs (eg. Young 
1980, Norman et al., 1985, Richstone k Tremaine 
1985), demonstrating that the observed structure of 
many galaxies is consistent with the presence of a 
MBH. At the same time, some authors have also 
shown that the observations may be accounted for 
by models of galaxies with no black holes (Duncan k 
Wheeler 1980, Binney k Mamon 1982; but see Mer- 
ritt 1987). 

When applied to this problem, most analytic tech- 
niques are restricted to spherical, or at best axisym- 
metric models of galaxies. Real galaxies are generally 
triaxial, and have small but measurable bulk rotation, 
which may strongly affect the influence a central MBH 
can have on the structure of the galaxy (Gerhard k 
Binney 1985, Pfenniger k de Zeeuw 1989). Using 
N-body realizations of galactic models, one can di- 
rectly examine the consequences of triaxiality, investi- 
gate instabilities (Merritt 1987, Palmer k Papaloizou 
1988), and analyze the orbital populations and ob- 
servational signatures of MBHs. In the limit of large 



N , simulations approach the intrinsic "graininess" of 
real galaxies, where the luminosity is supplied by a fi- 
nite number of effective point sources. Other authors 
have performed N-body simulations of the structure 
of galaxies with central MBHs, notably Norman, May 
k van Albada (1985) and Hasan k Norman (1990). 
As noted by Binney k Petit (1989), previous simu- 
lations (with N ~ 10 4 ) have been limited by poor 
resolution, spurious numerical relaxation or have em- 
ployed unrealistic distribution functions. 

Here we present a technique for simulating the adi- 
abatic growth of MBHs in galaxies, and discuss its ap- 
plication to spherical models. We compare our results 
with earlier theoretical results, in which MBHs were 
assumed to grow adiabatically in some background 
stellar distribution, and the final distribution func- 
tion was calculated assuming the action variables re- 
mained constant (Young 1980, Quinlan et al., 1994). 
The comparison serves to verify the code and analysis 
methods, validate the "adiabatic growth" approxima- 
tion for adding the BH to the galaxy, and check the 
final system for orbital stability in the presence of a 
BH. The code permits straightforward orbit classifi- 
cation for subsets of particles, and the analysis of the 
evolution of orbit families as the potential evolves. 

Using the Self-Consistent Field (SCF) method 
(Hernquist k Ostriker 1992), we can now follow large 
(N ~ 10 6 — 10 8 ) self-gravitating models for many 
dynamical times, allowing the central regions of the 
galaxies to be well-resolved, and making possible sta- 
tistically significant assertions about the spatial gradi- 
ents of observed properties of the models. The broad- 
band light profiles of elliptical galaxies are dominated 
by post-main sequence stars, constituting ~ 1% of 
the total number of stars. Our particle resolution 
is approaching the discreteness of the light tracing 
the potential; the actual potential is likely smoother, 
whether the potential is dominated by stars or dark 
matter. 

In this paper we examine how to generate and 
evolve large TV-body realizations of a family of distri- 
bution functions to which central MBHs are added. 
Using massively parallel processing (MPP) systems, 
we can integrate realistic calculations for ~ O(100) 
dynamical times, and analyze the intrinsic and pro- 
jected properties of the models as the MBH forms. 
We investigate the stability of the models, the lim- 
its of the assumption of adiabatic growth, and the 
true spatial resolution of our realizations as a func- 
tion of N and Mbh ■ We develop and verify schemes 
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for generating multi-mass realizations of our choice 
of distribution function, a multiple timestep scheme 
to allow fast integration despite the large range in in- 
trinsic timescales across our models, and parallelized 
analysis routines to calculate the projected moments 
of the distribution. 

2. Models and Methods 

2.1. The Self Consistent Field Method 

The simulation method we use is based on that 
described by Hernquist & Ostriker (1992). The po- 
tential of the galaxy, $(r, 9, <f>), is expanded in an or- 
thonormal set of basis functions. By choosing a suit- 
able biorthogonal basis set, the density, p(r,6,(f)), is 
represented by a similar expansion; namely 

p(r,6,<f>) = ^2 A nlm p nlm (r, 6, <f>) (2-1) 

nlm 
nlm 

where the A n \ m are the expansion coefficients for the 
basis chosen and 

V 2 $ n;m (r, 0, <f>) = 47r Pnlm (r, 9, <f>). (2-3) 

We choose as our standard zeroth order basis func- 
tion the Hernquist model (Hernquist 1990, Hernquist 
k Ostriker 1992), with units G = l,M = I, a = 1 
(note Quinlan et al., 1994 used a = 1/3), defined by 

^°oo(r) = / (2-4) 

27rr(l + ry 

$ooo(r) = T ^ : . (2-5) 

The particles determine the A n \ m through a nu- 
merical integration over the density, and move under 
the potential gradient derived from the expansion by 
equation (2-2), providing a completely self-consistent 
scheme for particle interaction. Each of the particles 
describing the density distribution may have different 
masses. The expansion coefficients can be saved dur- 
ing the time evolution and provide a compact snap- 
shot of the time evolution of the density distribution. 
Together with a sampling of the particle phase space, 
the A„i m provide a straightforward means for numer- 
ical classification of orbit families. 

Using this basis set we can well-represent a num- 
ber of standard potential-density pairs for spherical 



galaxies, including the truncated isothermal sphere 
and the family of 7-models (also known as ry-models) 
(Dehnen 1993, Tremaine et al., 1994). The Hernquist 
(or 7 = 1, r) = 2) model is our canonical example for a 
spherical galaxy containing no central MBH. All our 
results are compared with this default model, and the 
derived analytic results. 

The SCF algorithm is very parallelizable. An 
efficient parallel implementation of this technique 
has been implemented on some MPP architectures, 
specifically the CM-5 and T3D, with other implemen- 
tations under development (Hillis & Boghosian 1993, 
Hernquist et al., 1994). Using parallel SCF codes, we 
can integrate models with N ~ 10 6 — 10 8 for ;> 100 
dynamical times, on machines like the NCSA 512 pro- 
cessor CM-5. With 2 23 particles, typical of a full 
scale simulation, a file describing the complete par- 
ticle distribution (m;, a;;, ?/;, Zi,vxi,vyi,vzi,<!>i), re- 
quires 512Mb of disk space. Test models, such as 
we describe here, generally employ 512,000 particles. 
A typical test run, integrating a spherical model for 
~ 40 dynamical times, requires one or two hours using 
either 128 or 256 processors of the CM-5. 

2.2. The models 

We generate initial conditions from the distribu- 
tion function, f(E). The model is truncated at some 
radius, r t 1, and particle coordinates are chosen 
using an acceptance-rejection algorithm. For equal- 
mass models, in which m = 1/N, we have routines 
to generate realizations of various distribution func- 
tions, notable, 7 = 0, 7 = 1 (Hernquist) and 7 = 2 
(Jaffe) models (Hernquist 1990, Jaffe 1983). Given an 
initial realization, we grow a black hole at its center. 
We choose some black hole mass, Mbh , and a time to 
grow the black hole, tsn, and then introduce a mass, 
M(t), for t < t BH , 

Mm = M M ( 3 (J-) 3 -2(J-) 3 ),2. 6 , 

M{t) = M BH fort >t BH , (2-7) 
having an associated (softened) potential, 

*bh (r, t) = M(t)/y/r* + el H . (2-8) 

The particles move under the combined potential gra- 
dient, V($ + <&bh)- The choice of Ebh is driven by 
the spatial resolution in the center; higher N requires 
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smaller Ebh- For the equal-mass tests Ebh ranges 
from 2.5 x 10" 3 to 10" 2 . 

For a 7-model, the mass interior to some radius, 
r, is given by 

M(r)={j^)\ (2-9) 

In equal-mass models, the number of particles, inte- 
rior to r <C 1, N r , is then N r k, N x r 3_7 (Dehnen 
1993). For 7 = 1, and N = 10 6 , we have N r=0A « 
10 4 , and # r=0 .oi ~ 10 2 . So even with N ~ few x 10 6 , 
the model has no statistical resolution for r <J 10 -2 , 
just where we expect unambiguous signatures of a 
central MBH (Young 1980, Quinlan et al., 1994). 

To improve the resolution of our models we intro- 
duce a multi-mass scheme, where the mass of a par- 
ticle, m, is a continuous variable. We rewrite the dis- 
tribution function, f(E) = N(E , J) x m(E , J), where 
J is the particle angular momentum. Then we define 

* = < 2 - w > 

m(E,J) = m n r p for r p < r m , (2-11) 

where v t is the transverse velocity of the particle, M 
and a(= 1) are the total mass and scale radius, r m 
is some limiting radius (= 1 in practice), and m n is 
a normalizing scale factor. r p is an approximation to 
the particle's pericenter, usually accurate to within a 
factor of two. For r p > r m , m = const. A controls the 
range of masses used. 

In practice we calculate models for A = 0, 0.5, 0.75, 1.0. 
A = 0.5 provides a moderate mass range and A = 1.0 
provides a more extreme mass range. We suppress the 
mass variation outside r m so that the representation 
of the halo of the model remains tolerably smooth, 
and statistically robust. Figure 1 shows the mean 
particle mass as a function of r for A = 0.5, 1.0 in a 
7=1 model. The A = 1 model provides an order of 
magnitude increase in mass range over A = 0.5, and 
correspondingly larger numbers of particles at small 
radii. We tested multi-mass, 7=1, models for spuri- 
ous relaxation and mass segregation. No mass segre- 
gation was found, as might be expected by the nature 
of the force calculation, and any evolution was con- 
sistent with relaxation to virial equilibrium due to 
truncation of the initial conditions, and numerical re- 
laxation due to the finite number of particles. The 
relaxation was not significant for the large numbers 
of particles we employ in our simulations. 



With A = 1, we gain over two orders of magnitude 
in particle resolution near the center over equal-mass 
models. In order to make use of the improved spatial 
resolution, we are forced to a smaller Sbh 10 -3 , 
providing more than a factor of ten gain in spatial 
resolution. In the absence of a black hole, central ve- 
locities are ~ 1 in our units. With a central black hole, 
the velocity increases as r -1 / 2 down to the smoothing 
length, requiring a correspondingly smaller timestep 
for the integration. For a small Sbh ^ 10 -3 , and a 
large Mbh{^ 0.01), this forces a prohibitively large 
number of timesteps per dynamical time. To circum- 
vent this problem, we implemented a simple, efficient 
multiple timestep scheme. 

Ideally we want the particles near the black hole 
to move on smaller timesteps than particles at large 
radii. The particles requiring small timesteps con- 
stitute a near negligible fraction of the total mass, 
and the potential they are moving in is, in general, 
dominated by the black hole. The self-gravity of the 
rapidly moving particles is negligible. To retain par- 
allelism we want to avoid treating a subset of the par- 
ticles differently from the majority. We implemented 
a two level timestep scheme, whereby ^bh{i") is up- 
dated more often than $(r). We choose some control- 
ling timestep, dt, sufficient for a precise integration 
of the self-gravitating particles. Every dt, we recal- 
culate the expansion coefficients and the associated 
potential. In addition, we introduce a short timestep, 
dt s = dt/2 q . We update ^bh{i") for each particle ev- 
ery short timestep, and integrate the particle motion 
using a simple leapfrog integrator, under the updated 
black hole potential, holding $(r) fixed. For r <C a, 
where velocities are high, V$ ps V$bh, and at larger 
radii, the potential changes slowly compared to dt s . 
This scheme is very fast, q = 10 requires only a factor 
of 2 longer CPU time, while gaining a factor of 1024 
in time resolution for motion near the center, and a 
corresponding factor of 100 increase in spatial resolu- 
tion. Using this scheme, integrating 512,000 particles 
with no black hole present, and q = 10, energy is 
conserved to 8E/E ~ 10 -6 over 40 dynamical times, 
implying that the scheme is robust. 

2.3. Analyzing the data 

To analyze our simulations, we concentrate on two 
types of output. The intrinsic properties of the dis- 
tribution, such as volume density, dispersion, kurtosis 
and anisotropy, are useful for comparing with theoret- 
ical models. The projected properties of the model 
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must also be determined for direct comparisons with 
observations. We derive both the surface density and 
classical velocity moments, and the Gauss-Hermite 
moments of the projected velocity distribution, fol- 
lowing Gerhard (1993) and van der Marel & Franx 
(1993). 

For comparison to theory, we bin a model into an- 
nular zones, adjusting the width of the zones to con- 
tain equal numbers of particles. Summing over the 
mass weighted particle distribution in each zone, we 
derive the density and surface density, true and pro- 
jected dispersion, and the projected mean velocity, 
third and fourth moments of the velocity (skewness 
and kurtosis, k =< > / < vp 1 > 2 ), as well as the 
anisotropy, (3=1— < i> 2 > / < 2i> 2 >, where v r is the 
radial velocity. This choice of binning allows constant 
statistical sampling across the model, maximizing the 
signal for models with symmetry, and provides direct 
comparison with theory, specifically to the results of 
Young (1980) and Quinlan et al. (1994). Typically 
we use nj = 2000 particles per annular zone. 

In practice, observers do not fold their data into 
constant light annuli. For comparison with observa- 
tions, we sample our models with synthetic "slits" and 
apertures. The aperture sampling is simply done by 
considering all particles inside some projected radius 
Ro, and considering the projected properties of the 
distribution as Ro varies. 

For "slit" projection we consider the projected 
properties of a rectangular region, projected down the 
(arbitrary) z-axis of our model. The "slit" has some 
length, y s , and width, x s , divided into n s rectangular 
boxes along the y-axis. The slit is symmetric about 
the center along the y-axis but may be offset from 
the center along the x-axis by some value, x . For 
spherical models, slit analysis does not provide any 
additional information. The method was developed 
with consideration for future work where we will in- 
vestigate axisymmetric and triaxial models. 

For each box or annulus, we calculate a surface 
density, S(r), projected mean velocity, v z (r), pro- 
jected dispersion, cr(r), projected skewness and kurto- 
sis. In addition, we calculate the Gauss-Hermite mo- 
ments, Sj(r) and /i 8 (r) (Gerhard 1993, van der Marel 
& Franx 1993). Following Gerhard (1993), we define, 
wj(r) = (v zj - v z (r))/a(r), and 

I n 

s *( r ) = "A £ffiK-)e-* u ';, (2-12) 

II i=i 



where the Hi are the usual Hermite polynomials, and 
Pi = 1/V2 8_1 i! are normalizing constants. As noted 
by van der Marel et al. (1994), the observed velocity 
distribution is not fit for the true v z , a, but rather 
a "best fit" Gaussian profile is derived from the line 
profile. Hence they derive a "best fit" Gauss-Hermite 
fit, hi(r), defined as for the Si, but using "best fit" 
v' z (r) and cr'(r), such that hi(r) = = /12(f). We 
follow Heyl et al. (1994) and derive the /i 8 (r) moment 
coefficients iteratively from the s;(r) moments. Given 
v z ,<t,si,s 2 , define, v' zl (r) = v z (r), a[(r) = a(r), and 
solve for hi(r), then define 

= + h 1 x<r', (2-13) 

a' l+1 = *i + h 2 x*i, (2-14) 

and solve for hi,h 2 recursively until hi y2 (r) < 6. 
In this paper we choose 6 = 10 -6 , though the so- 
lution is not sensitive to the exact choice for 6, in 
general. About a dozen iterations are required for 
hi ^ to converge. The algorithm is easily paralleliz- 
able, and a version of the slit analysis has been im- 
plemented on the CM-5. Using logical parallel mask 
constructs on the phase space arrays, expensive sorts 
may be avoided and the data reduced rapidly. An 
N = 8,388,608 model is analyzed in less than 100 
seconds using 256 nodes of the CM-5. 

In the future, we also intend to generate synthetic 
line-profiles, using a blend of stellar line-profiles with 
appropriate offset, drawn from a library of model 
stellar lines. From such synthetic "observations" we 
can test how well analysis of actual observations can 
reproduce the underlying dynamics when convolved 
with seeing errors and observational noise. 

3. Results 

3.1. Tests of Parameters 

Preliminary tests of our code used the truncated 
isothermal sphere as the basic model. Comparison 
was made with the classic paper by Young (1980) and 
our results agreed both with those of Young and our 
separate re-analysis (Quinlan et al., 1994). 

In what follows our canonical test case is a 512,000 
particle equal-mass realization of a Hernquist model. 
All numerical results are compared with the basic re- 
sults derived from that model, and the correspond- 
ing analytic work by Quinlan et al. (1994). The 
model is truncated at r t = 300a. The mass enclosed, 
M{r) = r 2 /(l + rf (= 0.993 for r t = 300), is renor- 
malized to unity by uniformly rescaling the particle 
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masses. The resulting model is slightly sub-virial 
and settles into equilibrium in a few dynamical times. 
This transient evolution does not significantly affect 
the response to a growing central black hole, provided 
r-t is sufficiently large. For 7 = 2, the center of mass 
is significantly offset from the density maximum for 
N <j 10 6 , and the density cusp may be destroyed for 
r <J 10 -2 if the black hole is not correctly centered on 
the density cusp. The displacement of the center of 
mass is also an issue with multi-mass models, where 
the outer particles are more sparsely sampled at fixed 
N. 

We grow an MBH with mass Mbh = 0.01, smooth- 
ing length Ebh = 0.01, over tsn = 20 dynami- 
cal times (a dynamical time is defined naturally by 
t = a/v), using dt = 2.5 x 10 -3 . The model is allowed 
to settle for a further 20 dynamical times before the 
integration is terminated after 16,000 steps. For this 
reference spherical model, we use n = 16, / = m = 0. 
The properties of the initial and final distribution are 
shown in Figures 2 and 3. Figure 2 shows the charac- 
teristic Keplerian rise in dispersion due to the black 
hole. The cusp induced by the black hole is detectable 
for r <J 0.1. At this spatial resolution and Mbh there 
is no detectable anisotropy, which is consistent with 
the analytic predictions. Figure 3 shows k — 3 and 
/14. The numbers are in agreement with the results of 
Quinlan et al. (1994). Note that contrary to the case 
of the truncated isothermal sphere, both k — 3 and /14 
are flat in the presence of a black hole at this reso- 
lution, and rise in the absence of a black hole, again 
in agreement with analytic results. The calculated 
value of both k and /14 at the smallest R is limited by 
smoothing. 

3.1.1. Integration Parameters 

To test the robustness of our assumptions, we con- 
sider variations of the simulation parameters. 

We first try using larger timesteps, dt = 10 -2 , and 
find that the results agree with our canonical model, 
with an increase in [3 at the center at r = Ebh that is 
not statistically significant. Increasing dt by another 
factor of two leads to significant radial anisotropy due 
to integration errors as stars near the center are scat- 
tered by the black hole. 

We did a run with n = 32, to check that the 
model was adequately resolved with the canonical 
n = 16. There was no significant difference between 
the n = 16 and n = 32 models, implying that n = 16 



is sufficient. A smaller n would be adequate for most 
of the models considered here; the A n \ m suggest n k, 8 
would have sufficed for most of the runs, but for pur- 
poses of testing the code we chose to use n = 16. The 
variation of the kurtosis, k — 3, with r was a little 
smoother with n = 32, showing more clearly the peak 
in k near r = 3 x 10 -3 , but the difference was not 
statistically significant. 

It is somewhat surprising that the shape of the 
MBH induced cusp is not sensitive to n for r <C a. We 
reproduce analytic estimates for the cusp to a smaller 
spatial resolution than are sampled by the relatively 
low order radial expansion. The reason this approach 
works, is that at these scales the potential gradient 
is dominated by the black hole, not the self-gravity 
of the responding change in the stellar distribution, 
which is negligibly small in comparison. Our tests 
show that the code does correctly reproduce the true 
dynamics of the problem, down to r~ 2ebh- 

To test the stability of the model, we ran a sim- 
ulation with / = 6,dt = 10 -2 . There was no sig- 
nificant growth of low /, m modes; indeed for n = 
0,/ = m = 1 and / = l,m = the power in the 
modes decreased with time. There was large frac- 
tional variation in the coefficients for n = 0, / = 5, 6, 
m = 1, but the power was not significant due to nu- 
merical noise. There is some concern that Keplerian 
degeneracy of the fundamental modes of spherical sys- 
tems supports slow modes that violate the adiabatic 
approximation (Tremaine [personal communication], 
Weinberg 1994), but this does not appear to be a 
problem for our models. 

3.1.2. Model Parameters 

We varied tsn to test the validity of the adiabatic 
growth approximation, using tsn = 20, 10, 1, and 5 x 
10 -3 . When tsn is too small, a large anisotropy is in- 
duced in the distribution as central particles are scat- 
tered by the rapidly changing potential. A significant 
radial anisotropy is observed with tsn = 1, and for 
t BH = 5 x 10~ 3 (dt = 2.5 x 10~ 3 ), the anisotropy is 
maximal (P(r) — ► 1 for r — ► r t ). However, tsn = 10 
shows no significant spurious anisotropy compared 
to the run with tsn = 20, and we conclude that 
tBH ^ 10 is adequate for adiabatic growth approx- 
imation in spherical systems. 

This result was independently verified by looking 
at the change in radial action of orbits in different 
(fixed) spherical potentials, to which a time varying 
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Keplerian potential was added. The radial action was 
conserved to within ~ 1%, provided the time scale for 
the Keplerian potential to grow was ~ 5 — 10 radial 
periods (Quinlan et al., 1994). 

Finally, we grew smaller black holes, Mbh = 10 -3 , 
with t BH = 10, dt = 10 _2 ,2.5 x 10 -3 , and 2.5 x 
10 -4 ; for the last set, Ebh = 10 -3 . These models 
demonstrate the limitations of the equal-mass, single 
timestep models, even for N ~ 10 6 , and should be 
compared with the multi-mass models presented be- 
low. As Figures 4 and 5 show, for Mbh = 10 -3 there 
is little observable signature of the MBH at this res- 
olution. Even with 512,000 particles the signature of 
the MBH is barely noticeable. Decreasing Ebh allows 
the physics at small r to be explored further, but the 
finite number of particles leads to a loss of signal, and 
the required dt makes the simulation prohibitively ex- 
pensive. As Figure 5 shows, Ebh is critical to estimat- 
ing the kurtosis at the center of the model, while hi 
is a more robust estimator (cf. the behaviour of the 
short and long dashed lines). k(0) is very sensitive 
to a few high speed particles near r = 0, which are 
poorly sampled with a finite N , and not present for 
feasible Ebh in the equal-mass, single timestep inte- 
gration, hi is a more robust estimator, and shows 
the downturn at small radii expected from analytic 
calculations for this smoothing length. 

We also ran equal-mass models for 7 = 0,2, verify- 
ing the results in Quinlan et al. (1994) to within the 
resolution of the models. The 7 = model shows the 
expected steep density cusp. For the 7 = 2 model, our 
basis set cannot resolve the self-gravity of the central 
cusp well. Integrating a 7 = 2 model with no cen- 
tral black hole n = 16 and dt = 0.01, the final state 
deviated significantly from the initial model only for 
r <J 0.03 after 50 dynamical times. The MBH induced 
cusp for a N = 512, 000, 7 = 2, equal-mass model is 
limited less by the finite radial expansion than by Ebh 
and N . The steepening of the cusp due to the cen- 
tral MBH, seen in Quinlan et al., 1994, is observed to 
r ^ 2 x e B H, for M B h = 0.01, and e B h = 0.01. It 
is important that the black hole be centered on the 
density cusp and not the center of gravity of the dis- 
tribution. For 7 = 2, even with 512,000 particles, the 
density maximum may be offset from the center of 
mass by ;> Ebh, and the evolution of the central den- 
sity is to a flatter density profile, even with a central 
MBH added. The dispersion still shows an increase 
in the offset model, as it must for a locally virial dis- 
tribution. 



3.1.3. Multi-mass Models 

Using a A = 0.5 or 1.0 multi-mass model provides 
a dramatic improvement in spatial and mass resolu- 
tion. Figure 6 compares an N = 512, 000, A = 0.5 
model with a equal-mass model. Figure 7 compares 
the resolution of N = 512,000,£ Bff = 1 x 10 -3 , 
A = 1.0 and 0.5 models. The anisotropy predicted 
by analytic models (Goodman & Binney 1984, Quin- 
lan et al., 1994) is clearly evident in the multi-mass 
models for r ;> Ebh- For the A = 1.0 model, the spa- 
tial resolution is less than the softening length, the 
number of particles at r ^ Ebh is large and statisti- 
cally resolved on scales < Ebh- The spatial resolu- 
tion of the multi-mass model is a factor of 10 better 
than in the equal-mass case, and the dynamics can 
be followed to Ebh = 10 -3 with the same number 
of timesteps using the multistep scheme. In order to 
obtain statistically significant slit "observations" of 
the central Gauss-Hermite moments, the multi-mass 
scheme is critical, as the number of particles in each 
slit box is much smaller than for the annular projec- 
tion at fixed N and R. 

3.2. Velocity moments 

While we evaluate both s 8 (r) and /i 8 (r) (Gerhard 
1993, van der Marel & Franx 1993), in practice the 
calculated moments are equivalent for the cases con- 
sidered here. S4 and hi may differ if nj is small and 
the realization of the average line profile is poorly 
sampled, or when the line profile becomes highly non- 
Gaussian, and the higher moments (i > 6) are large. 
The hi moments have the virtue that hi = = /12 
by definition. By construction, /i 3 = for the models 
considered here, so all the information is contained in 
/14 (higher moments may be calculated, but at this 
resolution are too noisy to be of use). hi is anal- 
ogous to the kurtosis, but the exponential weighing 
suppresses the divergence of k that makes it a poor 
estimator (van der Marel & Franx 1993). As a first 
approximation, S4 ps hi, and k — 3 ps 8^6/14 for 
h 4 < 0.03. 

As can be seen in Figure 8, hi{r) by itself is not an 
estimator for central MBHs. With the addition of the 
MBH, the hi of the Hernquist model approaches that 
of the isothermal sphere without a MBH. The value of 
h 4 (R = 0) is limited by e B h = 10" 3 . For R ^ 2ebh 
the deviations of the velocity profile from a Gaus- 
sian are dominated by the intrinsic velocity profile of 
the underlying model. When averaging the velocity 
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profile over an aperture centered on the MBH, the 
smoothing and particle resolution at small R become 
very important. Both the smoothing of the potential, 
and the small number of particles leads to a deficiency 
in high velocity particles, which would lead to an in- 
crease in /14(C)) if properly included. 

Figure 9 explores the effect both of the smooth- 
ing and the particle resolution on the dispersion and 
/14 measured in a circular aperture centered on the 
galaxy, as a function of the aperture size, Ro- Ro 
may vary both because of improved instrument reso- 
lution for a particular galaxy, and from comparing 
galaxies at different distances. To mimic the true 
population of high velocity stars near the MBH we 
boost the particle velocities by a correction factor, 
v f z 1 — * v z x (\/ R 2 + £ B H I ' R) 1 / 2 , and compare with the 
uncorrected profile. As the particles were integrated 
in the smoothed potential their spatial distribution 
inside 2ebh is not correct in any case. In particular, 
the integration scheme undersamples the pericenters 
of particles on orbits near Ebh biasing the velocity 
profile to lower velocities. 

As can be seen in Figure 9, without the correction 
the best fit dispersion flattens at about 2ebh- With 
the velocity correction the best fit dispersion contin- 
ues to rise to the limit of the particle resolution. Over 
the range of apertures and models, the fit to a Gaus- 
sian line profile is surprisingly good, as can be seen 
from both /14 and the ratio of best fit dispersion to 
true dispersion. The bottom left panel illustrates why 
the kurtosis is a poor estimator of the velocity profile, 
the central value is sensitive to both the aperture size, 
the model resolution and the velocity error induced by 
smoothing the potential. For A = 0.5, and with no 
velocity correction, the central /14 declines with Ro- 
With the velocity correction the decline in /14 is at 
smaller Ro. With higher particle resolution, /14 rises 
at the smallest Ro well-sampled by the particles. As 
Figure 10 shows, this is entirely due to particles in 
the A = 1 model reaching smaller radii and broad- 
ening the velocity profile. For A = 0.5, even with 
the velocity correction, the number of particles inside 
Ebh is too small to raise /14 significantly. 

Figure 11 shows the slit analysis for a equal-mass 
N = 8, 388, 608 model and an N = 512, 000, A = 0.5 
model. The multi-mass model is critical to retain 
particle resolution into the center of the model, while 
a large total N is necessary to get statistically sig- 
nificant gradients for the moments of the distribution 
viewed through thin slits. At small y, the cusp in £ 



due to the presence of the black hole is clearly visible. 
The multi-mass model has relatively more particles 
at small y, and provides almost the same resolution 
in the central bin as the 8M model. The 512/; multi- 
mass model becomes very noisy at y ps 0.1, but tracks 
the larger model well at smaller y. This figure can be 
compared with Figures 4 and 5; the improved spatial 
resolution and higher particle number at small radii 
provides a clear signal of the surface density cusp and 
the change in dispersion and /14 due to the MBH. 

4. Conclusions 

The work presented here is primarily to verify the 
technique developed, by comparing the results with 
analytic and numerical calculations performed us- 
ing a completely different approach (Quinlan et al., 
1994). Within the range of variables where our mod- 
els are applicable, the results for intrinsic and pro- 
jected properties of the models agree with those of 
our previous paper. 

The models and analysis techniques are consistent 
with previous work, and other results in the literature 
(Young 1980, Goodman & Binney 1984, van der Marel 
1994a, b). The basic algorithms for realization of the 
models, integration and analysis are correct. 

Adiabatic growth is well approximated when black 
hole growth times of ;> 10 dynamical times are used, 
in accordance with independent semi-analytic esti- 
mates of the variation of the action of particle orbits 
in time varying potentials. With tsn long enough we 
can be confident that we are seeing the true response 
of the stellar distribution. We also find no radial or 
low /, m instabilities for these spherical models when 
a central black hole in the mass range 10 -3 — 10 _1 is 
added. 

For spherical models, it is clear that the central sur- 
face density profile provides a poor diagnostic of the 
presence of a central MBH, while the projected dis- 
persion provides a robust signal for the presence of an 
MBH, provided the observations have adequate spa- 
tial resolution. The refurbished Hubble Space Tele- 
scope is ideal for such observations of nearby galaxies. 

Gauss-Hermite moments are a promising statistic 
for constraining anisotropic and non-spherical mod- 
els with rising central dispersion, and may be used to 
discard extreme models producing rising central dis- 
persion in the absence of a large central dark mass ( eg. 
Duncan & Wheeler 1980). By themselves, the Gauss- 
Hermite moments do not provide a strong discrimi- 
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nant for a central MBH. They are primarily useful 
for constraining bulk rotation, anisotropy, departures 
from symmetry and the distribution of the underly- 
ing model given the projected properties. With the 
Gauss-Hermite moments, particularly hs and /14, we 
may preclude anisotropic models and triaxiality as the 
cause of rising projected central dispersion (Binney & 
Mamon 1982, Merritt 1987, van der Marel & Franx 
1993). Measurements of hi at larger radii can help dis- 
criminate between model distribution functions that 
produce similar surface density profiles. 

The 7=1 model is a good approximation to the 
surface density profile of elliptical galaxies over a wide 
range of radii (Hernquist 1990). Setiing the scale 
length a so that the effective radius, R e = 1.82a is 
similar to that for massive ellipticals (~ 3kpc), pro- 
vides a scale for the spatial resolution of our mod- 
els. With a softening length of 10 _3 a typical for our 
multi-mass models, we are resolving length scales of 
~ 1.5 pc. At a distance of 1 Mpc this corresponds 
to an angular resolution of ps 0.3". For comparison, 
M87 is at k, 16 Mpc, and the HST with 0.1" resolution 
can explore spatial scales of ~ 8pc, corresponding to 
length scales of ~ few x 10 -3 in our models. 

Physically, our model must break down in real 
galaxies for radii r <C 10 _3 a. Relaxation and stel- 
lar collisions become important at small radii in the 
presence of a cusp induced by a central black hole, 
so there is little point in exerting large computational 
effort to achieve much higher spatial resolution. It 
is still necessary to over-resolve the model. High ve- 
locity stars near the center contribute strongly to the 
velocity moments; and for triaxial models we expect 
particles on box orbits to explore the center of the 
model. 

In real galaxies, we would expect MBHs to form at 
the density maximum, but it is possible that the BH 
may wander away from the center of the galaxy by 
an amount comparable to the spatial resolution used 
here, and that the real density cusp formed is also 
flattened. For example, black hole growth may occur 
during a merger with a satellite galaxy, and the nu- 
cleus of the satellite galaxy may remain a distinct sub- 
system at r ;> a while the central MBH grows by gas 
accretion. In that event, reflex motion of the MBH 
relative to the orbit of the satellite nucleus can cause 
the MBH to wander about the density maximum by 
distances of order 10 pc on time scales of O(10 8 )y. 
Our models indicate that in this case the final density 
cusp should be flatter than if the central MBH stays 



fixed at the density maximum, possibly flatter than 
the underlying density profile at the break radius. 

Quinlan et al. (1994) found that the shape of the 
density cusp at the centers of galaxies is not a good di- 
agnostic of the presence of a central MBH. Steep cusps 
(eg. 7 = 2 models) may be present in the absence of a 
black hole, and the slope of a cusp induced by a cen- 
tral MBH may depend on the underlying stellar dis- 
tribution even at constant stellar density. The results 
from our simulations of "off-center" MBH growth re- 
inforces the conclusion that T,(R — ► 0) is a poor indi- 
cator for the presence of an MBH, and that galaxies 
may have shallow or no resolved cusps even with a 
central MBH present. 

It is perhaps surprising that the SCF method can 
be applied successfully to the problem of adiabatic 
growth of central masses in galaxies. The SCF pro- 
vides a relatively coarse spatial resolution when the 
numbers of expansion terms used is small. In fact, 
we find excellent agreement between our simulation 
results and those obtained independently by analytic 
calculations, since the self-gravity of the response of 
the stellar background is negligible in the inner re- 
gions where the dynamics are dominated by the cen- 
tral mass. Out tests support the earlier claims by, eg. 
Hozumi & Hernquist 1994 and Johnston & Hernquist 
1994, that the SCF technique will be a valuable tool 
for some problems in collisionless dynamics, and indi- 
cate that further theoretical studies, along lines sim- 
ilar to what we have presented here, are warranted. 

The methods described here can be applied to large 
N simulations of axisymmetric and triaxial models of 
galaxies. We hope to constrain self-consistent non- 
spherical models of galaxies containing central MBHs, 
and to survey the observable properties of families 
of different model galaxies containing MBHs. The 
method is also well suited for comparing particular, 
self-consistent realizations of theoretical models with 
the observed properties of individual galaxies. 

We are grateful to Greg Bryan for help with par- 
allelizing the SCF, and Jeremy Heyl for providing us 
with previously unpublished results. This work was 
supported in part by the National Center for Super- 
computing Applications, the Alfred P. Sloan Founda- 
tion, NASA Grant NAGW-2422 and the NSF under 
grants AST 90-18526, ASC 93-18185 and the Presi- 
dential Faculty Fellows Program. 
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Fig. 1. — The mean mass, m(r), vs radius, r, for two 
variable mass j = 1 models with A = 0.5 and 1.0. 
The plot shows a histogram of m for different radii, 
and curves giving cumulative particle number as a 
function of radius. The scale on the left y-axis refers 
to the histogramed data; the scale of the right y-axis 
refers to the connected points. 
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Fig. 2. — The properties of a 512,000 particle Hern- 
quist model with no black hole and with a Mbh = 
0.01 black hole, compared with the analytic results 
of Quinlan et al. (1994). Top left plot shows the 
volume density, with the cusp clearly seen at small r. 
The top right plot shows the surface density profile vs 
projected radius, R. The bottom left plot shows the 
dramatic rise in dispersion with the introduction of 
the black hole. The bottom right hand corner shows 
the anisotropy, [3. The plots are done using annular 
projection with nj = 2000 particles per annular zone 
unless otherwise stated. 
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Fig. 3. — The projected kurtosis (minus three) and 
/14 vs projected radius R for our canonical equal- 
mass 7 = 1 model with no central black hole (dot- 
ted lines) and with a Mbh = 0.01 central black hole 
(solid lines). The numbers shown in this figure were 
obtained using the annular projection discussed in the 
text. 
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Fig. 4. — The surface density and projected disper- 
sion for a 7 = 1 model, with a central black hole of 
mass Mbh = 10 -3 compared with the initial model, 
at three different resolutions. The short dashed line 
shows the model integrated with Sbh = 10 -2 , dt = 
10 -2 , the solid and long dashed lines show the same 
model with Sbh = 10~ 3 , dt = 2.5 x 10~ 4 , with differ- 
ent particle numbers per annulus. 
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Fig. 5. — As for Figures 4 but showing k — 3 and 
/14. Compare the behaviour of the short dashed curve 
and the long dashed curve for k — 3 and /14. For k — 3 
the smoothed, long rft integration approximates the 
no-MBH model. 



Fig. 6. — The surface density, projected dispersion, 
anisotropy and projected kurtosis for equal-mass and 
A = 0.5 multi-mass Hernquist models. The plots 
show the dramatic improvement in resolution with the 
multi-mass scheme. The slight offset of the (dashed) 
analytic curves is due to the finite width of the an- 
nular zones. The numerical calculations are centered 
on the mean radius of the zone, and the weighing of 
the averaged properties toward smaller radii within 
each zone leads to the offset observed compared to 
the analytically calculated properties. 
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Fig. 7. — As for Figure 6, but comparing A = 1.0 
and A = 0.5. At e B H = 10" 3 , the A = 1.0 over- 
resolves the central region. For non-spherical models, 
where the projection has to be done using "slits", the 
additional resolution is critical. 
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Fig. 8. — The fourth Gauss-Hermite moment for 
the isothermal sphere and Hernquist model with and 
without a central MBH. The upper panel shows /14 
for a equal-mass N = 512,000, 7 = 1 model and a 
N = 388, 660 equal-mass isothermal sphere. The ra- 
dius of the isothermal sphere has been scaled down 
by a factor of 10 to match the Hernquist model scale. 
The lower panel shows for a multi-mass, A = 0.5, 
N = 512,000, 7 = 1 model. 
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Fig. 9. — The effects of smoothing and particle num- 
ber on the line profile. The plots show variable aper- 
ture fits to the center of the galaxies, as a function 
of aperture radius, i?o- The innermost point is for 
the innermost 512 particles, the aperture varies in in- 
crements of 512 particles. The top left panel shows 
the best fit dispersion, a' p for four cases. The dot- 
ted line shows the A = 0.5 model, the long dashed 
line shows the A = 1.0 model. The short dash and 
solid lines show the same models, but with the veloc- 
ity boosted to its Keplerian value to compensate for 
smoothing. The top right panel shows the ratio of the 
best fit dispersion to the true dispersion. The bottom 
right panel shows the variation in /14. The bottom left 
panel shows k — 3 as a function of aperture size. 
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Fig. 10. — The line profile from the innermost pro- 
jected 10,000 particles, weighted by particle mass. 
The line profile strength, C, is shown with arbitrary 
normalization, with the £(A = 1) scaled by a factor 
of four to fit on the plot. 
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Fig. 11. — Slit projection of a equal-mass model with 
N = 2 23 = 8, 388, 608 particles, with no central black 
hole and a black hole of mass Mbh = 10 -3 , compared 
with the slit resolution of a A = 0.5, N = 512,000 
multi-mass model. The top left panel shows the num- 
ber of particles per bin, as a function of the slit posi- 
tion, y. The slit width was x s = 0.01a. The bottom 
left panel shows the resulting surface density (with 
arbitrary normalization). The top right panel shows 
the projected "best fit" dispersion, a' p as a function 
of y. The bottom right panel shows h^(y), showing 
the importance of large N to get reliable spatial gra- 
dients of projected moments with high resolution slit 
projections. 
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